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A novel approach to account for hard-body interactions in (overdamped) Brownian dynamics simulations is 
proposed for systems with non- vanishing force fields. The scheme exploits the analytically known transition 
probability for a Brownian particle on a one-dimensional half-line. The motion of a Brownian particle is 
decomposed into a component that is affected by hard-body interactions and into components that are 
unaffected. The hard-body interactions are incorporated by replacing the 'affected' component of motion by 
the evolution on a half-line. It is discussed under which circumstances this approach is justified. In particular, 
the algorithm is developed and formulated for systems with space-fixed obstacles and for systems comprising 
spherical particles. The validity and justification of the algorithm is investigated numerically by looking 
at exemplary model systems of soft matter, namely at colloids in flow fields and at protein interactions. 
Furthermore, a thorough discussion of properties of other heuristic algorithms is carried out. 



I. INTRODUCTION 

Many physical processes in soft matter systems, like 
colloids^— , polymers^ or biomolecular systems^, can be 
modeled in terms of Langevin equations where the in- 
fluence of the solvent is captured by friction and ran- 
dom forces^— The separation of the respective length 
and time scales of the solvent particles and the con- 
stituent particles of soft matter typically allows one to 
use the high friction limit where the velocities of the 
latter particles are strongly damped and inertia effects 
are thus negligible (overdamped Langevin equations). 
Another common idealization in the modeling of meso- 
scopic physical systems within the framework of Langevin 
equations consists in representing the extremely short- 
ranged and strong repulsive contact forces between par- 
ticles or between particles and space-fixed obstacles or 
walls by hard-body interactions. Important examples 
for physical processes which are investigated using these 
limits include the dynamical behavior of polymers and 
colloids in dilute as well as concentrated solutions, the 
self-assembly of biomolecules into biological structures 
like membranes or protein aggregates, the motion of 
molecules in biomolecular systems and the migration of 
macromolecules in microfluidic channels^. 

To unravel the properties of such physical processes, 
a thorough investigation of the particles' trajectories is 
often indispensable. To name just two examples, the 
knowledge of the typical trajectories is crucial to ana- 
lyze association pathways of complex-forming proteins or 
the migration mechanisms of biomolecules in microfluidic 
channels. As in most cases an analytical treatment of 
the Langevin equations of motion is not possible on such 
detailed, trajectory-wise level, efficient numerical inte- 
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gration schemes for solving these stochastic differential 
equations are essentiah n i 12 These schemes require a cer- 
tain degree of regularity of the involved force fields^ 
so that hard-body interactions (associated with reflec- 
tive boundaries^) can not be integrated directly due 
to their singular nature. Therefore, it is necessary to 
represent the effect of singular hard-body interactions 
in a tractable form. In order to avoid the numerically 
costly introduction of steep (and regular) auxiliary po- 
tentials, alternative potential-free methods seem partic- 
ularly promising with respect to computational efficiency. 
Such methods are based on detecting 'unphysical configu- 
rations', ie. 'collisions' or 'overlaps' between hard-body 
entities that arise from a regular numerical integration 
step, and on some rule that specifies how this unphys- 
ical configuration is to be corrected in order to create 
a physically valid configuration. Several heuristic meth- 
ods along these lines have been discussed in the liter- 
ature, ranging from Monte Carlo type schemes, where 
unphysical configurations are just rejected, to methods 
where overlapping particles with hard cores are displaced 
along their connection line. This is done by applying spe- 
cific rules that are based on geometric considerations or 
on the physics of elastic collisions*^— In recent years 
the emergence of event-driven molecular dynamics has 
inspired the development of trajectory-based heuristic 
'event-driven' Brownian dynamics schemes ! 17 ' 18 These 
methods often lack a thorough justification.iS, Apart from 
the rejection methodic and the event-driven scheme of 
Tao et alpi£ it is not obvious how the different methods 
can be applied to driven diffusive systems with additional 
non- vanishing force fields. 

In this work we propose a novel method to account for 
hard-body interactions in numerical Brownian dynamics 
simulations. The new approach reverts to the transition 
probability for the time evolution of a particle whenever 
a hard-body interaction has to be accounted for. This 
transition probability is determined by the Fokker-Planck 
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equation which provides an alternative way to describe 
the physics of stochastic processes^— Once this tran- 
sition probability is known, a physical trajectory of a 
particle can be directly generated. However, for many 
systems this transition probability is not known analyt- 
ically. In this paper we will argue and discuss how the 
analytical solution of the Fokker Planck equation for the 
one-dimensional half-line can be used to account for hard- 
body interactions in general, more-dimensional systems. 
The new algorithm has the advantage that it can be ap- 
plied to systems with spatially varying force-fields and 
can be used independently of the applied numerical inte- 
gration scheme. In addition, the effects of the introduced 
approximations due to applying the transition probabil- 
ity for a system on a half-line can be critically assessed 
so that these approximations can be justified and numer- 
ically controlled. 

The article is organized as follows. In Sec. [H]the ba- 
sics of Brownian dynamics simulations are recapitulated. 
Then the problem of accounting for hard-body interac- 
tions in Brownian dynamics is discussed and the basic 
idea of our novel approach to incorporate them is out- 
lined. In this context other proposed heuristic methods 
are discussed in Sec. IIIII Section llVl is then devoted to the 
development of the new algorithm for space-fixed hard 
walls. The algorithm is next generalized to systems with 
hard spheres in Sec. [Vj The validity, quality and effi- 
ciency of the methods are scrutinized in Sec. I VII by nu- 
merically investigating different exemplary physical situ- 
ations with non- vanishing force fields. A brief account of 
some of the results of this article was already published 
elsewhere i2£ 



II. HARD-BODY INTERACTIONS IN BROWNIAN 
DYNAMICS 

A. Basics of Brownian Dynamics simulations 

Consider a Brownian particle in a solution and in 
a force field F = F(r, t) at temperature T. Apart 
from the deterministic force F due to external poten- 
tials, non-conservative forces or inter-particle interac- 
tions, the Brownian particle experiences the influence 
of the medium. One influence is related to the friction 
which is determined by Stokes' friction force — rjr with 
friction coefficient r\ and particle velocity r. In addition, 
the medium influences the particle by random collisions 
through its constituents. This thermal heat bath is repre- 
sented by a random force £ with the statistical properties 
of Gaussian random variables: (£ a (t)%P (s)) = S a ^6(t — s) 
with a and j3 denoting Cartesian components. In the 
strong friction limit, for which inertia effects are negligi- 
ble, the evolution of the particle is then governed by the 
Langevin equation 

f(t) = 5*1 + VW£(t), (i) 



where the strength of the fluctuating force is related to 
the diffusion constant D = k B T/rj (overdamped Langevin 
dynamics, Brownian dynamics) ; 3 ' 7 ' 12 Note that the as- 
sumption of the strong friction limit is well satisfied for 
biological, colloidal or polymeric systems4£ 

The Langevin evolution equation can be used to obtain 
the trajectory of the particle whose physical contents can 
then be investigated. In this context the Langevin equa- 
tion is best formulated in terms of differentials, 

dr = vdt + V2DdW. (2) 

Here the drift term is represented in terms of a veloc- 
ity v = F/rj and the stochastic term is related to a 
Wiener process W. The Wiener process W describes free 
Brownian diffusion and is mathematically defined to be a 
Gaussian stochastic process with mean (W a (t)) = and 
covariance {W a {t)W f3 {s)) = 5 afi min(i, s)£21i22 In ad- 
dition, the ltd Lemma (dW a (t)dW p (t)) = dt5 afi holds 
true for its differential. Note that the incremental Wiener 
process dW{t) cannot be written in the form w(t)dt due 
to the non-differentiability of the trajectory of a freely 
diffusing particle ^ 21 ' 22 

The numerical integration of @ will generate a dis- 
crete approximation of the trajectory. To this end a dis- 
crete 'incremental' displacement corresponding to © is 
generated for a finite time step from the current phys- 
ical position of the particle. The most prominent inte- 
gration scheme for such Brownian dynamics simulations 
is the Euler algorithm (see Sec. IIII[) . Apart from single 
step algorithms like the Euler method, multistep inte- 
gration schemes such as the Heun algorithm have been 
developed that additionally incorporate the information 
from configurations for intermediate steps (more details 
about the numerical integration of stochastic differen- 
tial equations, in particular on necessary conditions on 
the regularity of the involved forces, can be found in the 
literatureSJii 2 .^). 

One crucial point in carrying out a numerical integra- 
tion of © is the choice of an appropriate size of the time 
step between two successive times of the discrete trajec- 
tory. For a non-zero and spatially varying force field F, 
the time step has to be small enough so that the force 
is virtually constant on the scale that is associated with 
the typically proposed spatial displacement. For systems 
whose time evolution is usually modelled within a Brown- 
ian dynamics approach, however, the corresponding force 
fields are in general relatively smooth. This allows the 
use of fairly large integration time steps. 

B. Algorithm for accounting hard-body interactions 

For systems with hard-body interactions the direct nu- 
merical integration of the unmodified Langevin equations 
over a given time step can result in unphysical configu- 
rations where the particles are displaced into the wall, 
for example, or where hard-core particles penetrate each 



3 



other. Such configurations, which are obtained from the 
direct integration of the Langevin equation but are not 
consistent with hard-body interactions will be called 'un- 
physical' in the following. As the hard-body interactions 
are singular in nature they cannot directly be incorpo- 
rated into the integration scheme. One possibility is given 
by a representation in terms of regular potentials. This 
approach, however, has the drawback, that the auxiliary 
potentials have to be rather steep and therefore the ac- 
ceptable time step for the numerical integration has to 
be very small rendering the simulation extremely time- 
consuming. Thus other more effective potential-free ways 
to account for hard-body interactions are desirable^ 

The method we propose here is based on the approach 
to consider the physics of a diffusing particle in a force 
field F by directly looking at the transition probability 
p(r,t;f ,t ) to find the particle at position r at time t 
provided it was at position fo at time to < t. This (nor- 
malized) probability is governed by the Fokker-Planck 
equation^— 



dt 



p(f,t\f Q ,tQ) = DA? 



- F (r,t)\ , 

V? ] p(r,t;r ,t ). 



(3) 

To solve this partial differential equation boundary con- 
ditions for the physical domain have to be provided. The 
presence of hard walls at a boundary with normal vector 
ft corresponds to reflective boundary conditions for which 
the normal projection ft ■ j of the probability flux 



j(f, f; r , to) = - ( DV P - J p(f, t; f , to) (4) 



vanishes at the boundary. 

In general, a realization for the displacement f — fo of a 
Brownian particle after elapsed time t — to can be directly 
generated from the transition probability p(f, t; fo, to). 
Suppose now that this probability is known in the vicin- 
ity of the hard surface and that the numerical integration 
using the unmodified Langevin equation (without hard 
walls) has led to an unphysical position r for the particle 
that was initially at fo. Instead of applying a modified 
integration step to the Langevin equation, a new physi- 
cal position r * for the elapsed time t is directly generated 
from that part of the transition probability p(r*, t; fo, 0) 
that represents the particle's encounter with the wall. 
This recourse to the transition probability whenever a 
hard-body interaction has to be accounted for forms the 
heart of our proposed algorithm. In the following we will 
demonstrate that the transition probability p(f, t; fo, to) 
for spherical particles can indeed be decomposed in a 
mathematically rigorous way into a 'free' part (without 
hard-body interactions) and a part that stems from the 
reflective hard boundary. This decomposition is based 
on the solution of the Fokker-Planck equation on a one- 
dimensional half-line with a reflecting boundary at the 
origin. 



C. Exact one-dimensional algorithm 

For the particular example of a one-dimensional sys- 
tem, where the position q is restricted to the half-line, 
[0,oo [ and is subject to the constant force / the Fokker- 
Planck equation reads 







d 2 fd 



dt P ^ q dq 2 ^ J] dq'' 



(- r >) 



The transition probability can be calculated analytically 
for a reflective boundary at q = . 25 i 26 Using the notation 
V := f/rj and setting to = without loss of generality 
the solution for the transition probability p(q,t\qo) = 
p(q, t; g , t = 0) is given by 

p(q,t;qo) =Pi(q,t;qo) +p 2 (q,t;q ) + p 3 (q,t;q Q ) (6) 
with the contributions 



pi{q,t;qo) 



P2(q,t;q ) 



exp 



(g-go- vtf 
4D q t 



, (7a) 



y/^Dql 



and 



P3(q,t;qo) 



2D, 



■ exp 



vq 
~D~„ 



erfc 



q + q + vt 



y/4D~i 



(.7c) 



Each pi satisfies the Fokker-Planck equation ([5]) indi- 
vidually, their sum then satisfies the initial condition 
p{q,t = 0;qo) = S(q — q®), the reflective boundary con- 
dition j(q,t;q )\ q = = -{D q -§^p ~ vp)\ q=a = 0, and is 
normalized on the half-line [0,oo[. The different parts 
have the following intuitive interpretation. The first term 
Pi is the probability distribution that corresponds to un- 
restricted diffusion in a constant force field and repre- 
sents the contributions from particles that do not en- 
counter the wall at q = during the propagation over 
the time t. The propagation of such particles is al- 
ready correctly captured by the 'free' Langevin equa- 
tion dq = vdt + y / 2D q dW corresponding to (|5|) without 
boundary at q = 0. The second and third terms contain 
the contributions from particles that encounter the wall 
during the propagation time t. These contributions p 2 
and P3 are obviously not correctly contained in the asso- 
ciated Langevin equation without further consideration 
of the reflective boundary condition at q = 0. 



1. The algorithm 

It is now obvious how physical positions that are con- 
sistent with a reflective boundary at q = can be gener- 
ated: 

1. Starting at q > 0; suggest a new position q of 
the Brownian particle (with diffusion constant D q ) 
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after time step dt by some given numerical inte- one can write pn = (p 2 +P3)/w. The random variable A 

gration scheme for the Langevin equation (without associated with the algorithm has the distribution func- 

further consideration of the hard wall). tion 

2 a. If the proposed q is in the physical domain, i.e. +OG oo 

q > 0, it is accepted as the new particle position. , \ f , f , , N , N /ir .x 

Paigl?) = / dxi / dx 2 p G {x x )p^{x 2 ) (12) 

2 b. If the proposed position is unphysical, i.e. q < -°° o 

0, a new position is generated directly from the xS(q— [x\9(xi) + x 2 0(— Xi)]) (13) 

transition probability (see Appendix [XJ -foe oo 

„ = / /d gaPl (x 1 ) ^+^ (14) 

Jo dg +p 3 (?)) x<y(g- [0:^(0:1) +^(-3!!)]). (15) 



Obviously, the end positions which are directly accepted 
from the integration step of the Langevin equation make 
up the contribution pi in the exact solution ([6]) of the 
transition probability. The distribution pn represents the 
(normalized) part p 2 + P3 in the analytical solution ^ 
capturing the physics of particles that encounter the wall 
at least once. See Appendix [A] on how to generate q from 
(|8|) in practice. 



2. Discussion 

The algorithm is constructed such that a stochastic 
particle trajectory is generated from two contributions: 
Parts where the particle does not encounter the wall at 
all are captured by a standard integration step, whereas 
the parts of the trajectory where wall-collisions occur are 
constructed from pn . Let us look at this point from a dif- 
ferent angle. Mathematically speaking the algorithm uses 
two random variables A G and Ah to construct a valid end 
position X. The Gaussian variable A G describes unre- 
stricted free diffusion, can thus take on all real numbers 
R and is distributed according to the Gaussian distribu- 
tion p\ (compare (|7a|)). which is normalized on R. Ah 
can take on only positive real numbers R>o and is dis- 
tributed according to (the normalized) pn of relation |[8} . 
The final random variable generated from the algorithm 
is of the form 

A = A G 0(A G ) + A H 0(-A G ), (9) 

where 9{X) is the Heaviside step function. Defining the 
collision probability 

00 

W = J dq{p 2 (q) + p 3 {q)) (10) 


00 

= Jdq(l- Pl (q)) = J dgpi(ff) 

-00 

= lcrfcf g0 ±^ (11) 
2 \^dt) 



Due to the Heaviside functions in the argument of the 
delta function the latter can be rewritten as 5(q — 
[x 1 8(x 1 )+x 2 6(-x 1 )}) = e{x 1 )5{q-x 1 )+e{-xi)5{q-x 2 ). 
The resulting two integrals finally add up to Paig(?) = 
Pi(q) +^2(9) +P3(q) and hence the algorithm indeed gen- 
erates random variables which are distributed according 
to the Smolouchowski solution (0- 

By construction our proposed procedure for a simula- 
tion in a system with constant force will generate a tra- 
jectory with the correct statistics independent of the size 
of the used time step dt. For a slowly varying force the 
time step can be chosen small enough so that the force 
stays virtually constant for the typical displacement of 
the particle. For each time step the simulation is then 
carried out as above with the presently acting force (see 
Appendix IB] for an illustration). For such an appropri- 
ately adjusted time step the Brownian dynamics simula- 
tion will generate a trajectory that exhibits the correct 
physical properties in a controllable accuracy. 

For the outlined approach to work for a general phys- 
ical system one has to solve the Fokker-Planck equation 
for the transition probability in the vicinity of a hard 
surface analytically. This is only possible in closed form 
for a very limited number of problems. In the follow- 
ing sections we will show in details that the solution for 
the one-dimensional half-line can be used to construct 
an approximate transition probability with which hard- 
body interactions can be accounted for in many situa- 
tions. The approximation is thereby controlled by the 
size of the time step. For small time steps only small 
displacements show up so that a generally shaped hard 
surface appears to be virtually flat. For the normal direc- 
tion the systems is then basically restricted to a half-line 
whereas the motion in the components parallel to the 
surface is unrestricted. This parallel motion is decoupled 
from the normal motion and can thus be treated within 
the usual approach whereas the one-dimensional normal 
motion is dealt with by the one-dimensional approach 
from above. Details will be worked out below. We note 
already at this point, however, that the involved adaption 
of the time step for spatially structured hard boundaries 
works along the same lines as the corresponding choice 
to capture the effect of spatially varying force fields. 
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III. HEURISTIC REPRESENTATIONS OF HARD WALLS 

Before extending the outlined approach to more com- 
plex physical systems we discuss the quality and valid- 
ity of other proposed heuristic methods to account for 
hard-body interactions in a potential-free manner. To 
this end we consider a Brownian particle on the one- 
dimensional half-line in a constant force field for which 
the analytical solution © is at hand. This special sys- 
tem will enable us to precisely identify the shortcomings 
of those schemes. In the following we will only compare 
methods that can be directly extended to systems with 
non- vanishing forces. 



A. Rejection scheme 

A first possible scheme could consist in discarding all 
displacements into hard walls as unphysical and repeat- 
ing the integration step until a valid displacement is ob- 
tained. However, in doing so the influence of the contri- 
butions P2 and P3 is not captured even approximately, as 
this scheme tries to reconstructed the whole transition 
probability only from (the renormalized) p\ . This results 
in a very poor convergence in the limit of decreasing inte- 
gration time step dt as can be seen in Fig. [1] (left panel). 
In particular, one would expect that systematic devia- 
tions show up only in close proximity of the wall. The 
mentioned renormalizing of p\, however, leads to devia- 
tions even far away from the wall (compare Fig. [TJ left 
panel). 

The Monte Carlo inspired rejection scheme models the 
effect of hard body interactions by dismissing a proposed 
propagation that leads to an unphysical configuration, 
but still advancing time.— The part p\ in the probability 
distribution (JSJ) is correctly reproduced for a time step 
dt, but the contributions from pi and p^, are replaced by 
a delta function situated at the initial position qo with a 
weight given by the collision probability w (compare (|10| . 
see also Fig. Q] (right panel) and Fig. [5]). The one-step 
distribution function of the rejection scheme thus reads 

Pre(g) =Pi(q) + ^erfc ^^= \ S(q-q ). (16) 

In this way the relative weight of the contributions p\ and 
Vi + P3 is correctly captured, although the influence of 
encounter events is not 'explicitly' worked in. In addition, 
the reflective boundary condition at q = is violated. 
Instead of a vanishing probability flux one has 



cxp 



Jrc\q=Q — 



( ( g0 +t'dt) 2 \ 
V iD q dt ) 



qo — vdt 



^AnDqdt 2dt 



(17) 



In practise a finite time step will give acceptable re- 
sults for large initial separations from the wall where the 
typical displacement will not lead to an encounter. Close 



to the wall where yj2D q dt/ qo is considerably larger than 
zero, however, systematic deviations show up. For illus- 
tration we give (the asymptotic series of) the mean value 

(<?)rc of (mi) 



(q)re = (q) - 



^AnDqdt 



e 2D i2D q dt + 0(dt z )) (18) 



where (g)is the mean of the exact solution ©. (In Refs. 
[3 and M the convergence of the rejection method in the 
limit of vanishing dt was discussed for systems without 
deterministic forces). Figure [1] (right panel) illustrates 
the influence of the size of the time step on the conver- 
gence behavior. The chosen parameters correspond to a 
situation where the probability to encounter the wall is 
already moderately small. Even then rather small time 
steps have to be used to reduce the systematic deviations 
for end positions close to the wall. 

We note that a modification of the rejection scheme 
was proposed in Ref. [27] for many-particle systems. 
There the amount for which time is advanced during 
one step is related to the fraction of proposed displace- 
ments without unphysical hard-core overlaps compared 
to displacements with overlaps (rejection with 'intrinsic' 
clock). For just one particle, this scheme corresponds to 
the rejection method without advancing time, which has 
been demonstrated above (see also Fig. [TJ left panel) to 
exhibit poor convergence. 



B. Event-driven methods 

The event-driven scheme of Tao et al.— attempts to 
capture the influence of encounter events by using an 
auxiliary path for the propagating particle. For an in- 
tegration step where an overlap is detected the singular 
trajectory of a Brownian particle is replaced by a regular 
path whose parameterization is obtained from the used 
integration scheme. To keep the discussion clear we con- 
sider the Euler algorithm where the new particle position 
q(t + dt) is given by 

q(t + dt) = q(t) + v{t)dt + yj2D q dtG (19) 

with G being a random Gaussian variable of zero mean 
and unit variance. Suppose the particle starts from a po- 
sition q(t) = qo close to the wall and its end position is 
in the wall. The regular auxiliary path is then parame- 
terized by q(t + Xdt) = q + vXdt + y/2D q dt\f\G with 
A € [0, 1] and fixed realization G. From this path a 'col- 
lision time' Ao with the wall can be calculated (note that 
Ao itself is a random variable) . Then the particle is prop- 
agated along the path to the collision point corresponding 
to an elapsed time Aod£. Starting at the wall the particle 
is then propagated the remaining time dr = dt — Xodt 
with a new realization of the Gaussian random variable 
such that the initial part of the path leads away from the 
surface and a new physical end position is obtained. This 
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FIG. 1. Transition probability of a particle on a half-line for 
the parameters v — —1.0, D q = 1.0, initial position go = 0.5 
(corresponding to an encounter probability w, equation (|10|) . 
of approximately 8%) and elapsed time t = 0.05:— The ana- 
lytical curve ([6]) is displayed by the solid line and is compared 
to results from numerical simulations based on the Euler al- 
gorithm using the rejection method without advancing time 
(left) and updating time (right). The (red) dotted curve is ob- 
tained with dt = 0.05. For the rejection scheme with updating 
time a delta-like peak at the initial position emerges repre- 
senting the probability to encounter the wall.— The (blue) 
dashed and the dotted-dashed curves were obtained using 
dt = 0.005 = t/10 and dt = 0.0005 = i/100, respectively. 
(For comparison the light (orange) line in the inset shows 
results from the event-driven method discussed below with 
dt = t = 0.05.) 



approach obviously captures the influence of encounter 
events, however, the contributions from pi and ps are 
not completely collected. Suppose, for instance, that the 
deterministic force is pointing away from the wall, i.e. 
v > 0. For a negative Gaussian variable G the path for 
the propagation for the remaining time dr starting di- 
rectly at the wall has an initial part which will always 
lead into the wall although it might produce a valid end- 
position. Such a path has to be rejected. A physical end 
position q £ [0, oo[ is therefore only generated if the sign 
of the Gaussian random variable G is positive. The phys- 
ical end-position is therefore systematically shifted away 
from the wall by an amount of 0(vdr). To avoid this 
systematic shift one can relax the conditions on the ran- 
dom Gaussian variable and allow all paths that lead to 
valid end-positions. We use this latter variant for the dis- 
cussions throughout the present paper. In addition, the 
event-driven approach suffers from the serious drawback 
that the sampling of a probability distribution directly 
from trajectories never can generate the negative contri- 
butions of P3 for forces v > 0. Hence, a procedure that 
tries to incorporate the presence of a hard wall into a 
scheme that generates (discretized and regular approxi- 
mations of the) trajectories will not correctly capture the 
contributions of p%. 

Figure [5] displays the analytical solution © of the 



Fokker-Plank equation © together with transition prob- 
abilities obtained from the event-driven scheme for differ- 
ent initial separations go from the wall. As expected, the 
systematic deviations are reduced for increasing go for 
which the encounter probability is diminishing as well. 

In principle the one-step distribution function for the 
event-driven scheme can be constructed analogously to 
© from X = X G 9(X G ) + X K 6(-X G ) where X K rep- 
resents the random variable generated in case of an en- 
counter. This variable is distributed according to 



exQ(- {x ~ vdT} ' ) 



p K (x;dr) 



v / ^^( 1 + erf( v /^))' 



x > 0. 



(20) 

However, the time dr depends on the value of X G < of 
the first integration attempt and is determined by 



X G = q + vdt + y/2D q dtG, 



(21) 



= qo + v{dt-dT) + ^j2D q {dt-dT)G, (22) 

from which G can be eliminated and the physical dr G 
[0,dt] be computed. The resulting distribution is thus 



Pev(q) =Pi(q) + / dxp K (q;dT(x))p G (x) . (23) 



This expression, however, cannot be evaluated in closed 
form. Recall that the event-driven scheme produces a 
distribution of remaining propagation times dr G [0,dt]. 
Due to relations (|2ip and (|22|) this can be expressed in 
terms of unphysical x G] — oo, 0] as is done in the integral 
in ((231). 

We already noted that the event-driven scheme in- 
cludes effects due to encounter events with the wall. In 
particular, all events that encounter the wall precisely 
once are correctly captured whereas those with several 
encounters are approximated by trajectories with one en- 
counter only. This observation readily suggests a refine- 
ment of the event-driven scheme. Once the particle is 
propagated to the first collision point the event-driven 
scheme is applied to the remaining time interval so that 
it is either propagated directly into the valid region or 
again a further collision point and collision time have to 
be determined. This approach is then iteratively applied 
till the particle is propagated for the total time step dt. 
From numerical investigations (not shown in this arti- 
cle) we conclude that this iterative event-driven scheme 
indeed leads to correct transition probabilities for forces 
into the wall. For forces away from the wall, however, the 
procedure is not correct as the negative contribution from 
P3 cannot be captured in this trajectory-based method. 

We note at this point that Lamm and Schulten2£ de- 
veloped also a scheme which successively takes overlap 
events into account. If an overlap is detected a further 
displacement away from the wall is determined with a re- 
maining propagation time that is sampled with the help 
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ing particle on a half-line subject to a constant force: 



FIG. 2. Transition probability of a particle on the half-line for 
v = 1.0, D q = 1.0 and elapsed time t = 0.05.— Three differ- 
ent initial positions qo are shown, namely, from left to right, 
qo = 0.05 (corresponding to a probability w of approximately 
38% to encounter the wall), qo = 0.35 (10%) and go = 0.6 
(2%). The solid curves depict the analytical solutions 
The numerical results are obtained by one integration step of 
the Euler scheme. The (red) dashed curves correspond to the 
event-driven method. For comparison the (blue) dotted line 
represents the results from the rejection method. In these 
curves a delta-like peak shows up at the initial position qo 
representing the particles that encounter the wall, and are 
thus not displaced in the rejection method (indicated by a 
spike) . The insets show the deviation of the results from the 
event-driven scheme close to the wall, each with the same 
magnification. 



of an approximation of the collision time distribution. If 
this still leads to an overlap a third step is generated 
according to (|7c[) approximated by a Boltzmann distri- 
bution. For details see Ref. |30L 



fc„(dt;g ) 



i\dt 



dtt n p(qo + (,dt;q ). 



(24) 



We remark that relation (|24|) is used to estimate 
the drift term and the diffusion coefficient in Fokker- 
Planck or Langevin equations, respectively, from empir- 
ical data i 31 ' 32 For go > one has limdi-s-o fci(dt; go) = v, 
lim dt ^ k 2 {dt; q ) = D q and \im dt ^ k n (dt; q ) = for 
all n > 2, i. e. Gaussian statistics is recovered in the 
limit of vanishing time steps dt. For finite dt, the second 
and third coefficients are shown in Fig. [3] as a function 
of the initial separation from the wall for the different 
integration schemes. The moments exhibit a strong de- 
pendence on qo close to the wall and show a non-Gaussian 
behavior. For example, the third moment 'vanishes' only 
for considerable distances from the wall. In general, for 
increasing distance go the coefficients approach the Gaus- 
sian values for unrestricted diffusion in a constant force 
field ('bulk-like' situation), namely D q + ^v 2 dt for k 2 and 



D q vdt 



0(dt 2 ) for fc 3 . 




C. Supplementary remarks 

Due to the presence of the hard wall (reflecting bound- 
ary) at the origin, the statistical properties of the par- 
ticle displacements for finite propagation times dt devi- 
ate from a Gaussian (compare the Smoluchowski solution 
©)• For a diffusing particle close to a wall, for instance, 
the first central moment (q — qo) will not be equal to vdt 
and thus will not vanish even for absent deterministic 
drift. Only in the limit ^2D q dt/q (and g > 0) 
the 'Gaussian' statistics is restored, i. e. for small prop- 
agation time and/or sufficiently far away from the wall, 
so that its influence is practically not noticeable. For 
analyzing how well the different integrations schemes re- 
produce these deviations from a Gaussian behavior close 
to the wall, we quantify them by looking at the nth ex- 
pansion coefficient of the Kramers- Moyal series^ - — that is 
associated with the transition probability © of a diffus- 



FIG. 3. Second and third ('finite-time') coefficient fe(dt;go) 
(lower left panel) and k3(dt;qo) (lower right panel) of the 
Kramers-Moyal expansion as defined in Q24p for a diffusing 
particle on a half-line with v — 1.0, D q — 1.0 and elapsed time 
dt = 0.05. 28 The solid curves depict the k n from the analytical 
solution The results are obtained using one integration 
step of the Euler scheme. The (red) dashed curves correspond 
to the event-driven method (using (|23|) and numerical integra- 
tion), the (blue) dotted line represents the rejection method 
(with advancing time). The event-driven data is only slightly 
off the exact curve (compare insets). The upper graphs on 
both sides show the probability w for an encounter with the 
wall, see relation (|10p . 

From Fig. [3] we see that in particular close to the wall 
the rejection scheme shows significant deviations from the 
exact solution. Its convergence properties do not only de- 
pend on the choice of dt with respect to the force field, 
but also on how close a trajectory approaches the wall. 
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Figure 2] illustrates this mutual dependence of a proper 
dt and the distance by analyzing how good the statisti- 
cal properties of the stochastic evolution are captured. In 
contrast, the choice of dt for our algorithm is not influ- 
enced by the distance from the wall and is solely limited 
by the variations in the force field. 
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FIG. 4. Second and third ('finite-time') coefficient fc2(dt;go) 
(left hand side) and ks(dt;qo) (right hand side) of the 
Kramers- Moyal expansion as defined in (|24[) for the differ- 
ent algorithms as a function of the time step dt (v — 1.0, 
D q — 1.0) to illustrate their convergence.™ Suppose the par- 
ticle to be propagated by dt is located at go = 0.5. The 
coefficients and hj, are displayed in the upper panels. The 
black curves represent the exact expressions (by construc- 
tion identical to our algorithm), the (red) dashed curves the 
event-driven scheme and the (blue) dotted curves the rejec- 
tion method. (The light grey curves show the asymptotic 
expressions for dt — > 0, 'bulk-like' situation, see main text). 
For the choice dt = 0.01 (vertical dashed lines) ki and kz are 
well reproduced by the rejection scheme. Suppose now that 
the particle is propagated to q = 0.35. The coefficients for the 
next time step starting then at qo = 0.35 are displayed in the 
lower panels. Now, the rejection scheme shows considerable 
deviation in k% and &3 for the originally adjusted dt = 0.01. 



that the time step is already chosen appropriately small 
to account for variations in the force fields. 



A. Description of algorithm 

Suppose that the Brownian particle is currently at po- 
sition r close to the wall and experiences the deterministic 
force F. Then, the standard integration scheme, which 
'suggests' a displacement dr, is modified in the following 
way to yield a physically valid displacement dr*, which 
accounts for the hard wall: 

1. Calculate a 'preliminary' displacement dr using the 
standard integrator for the 'free' Langevin equation 
without wall. 

2. Check whether or not the so suggested new position 
r + dr yields a physically valid configuration (i.e. 
whether or not the particle, if at position r + dr, 
overlaps with the hard wall) . 

a. If this configuration is valid, the suggested dis- 
placement is accepted as the new one, dr * = 
dr. 

b. If this configuration is invalid, the new, valid 
displacement is calculated according to 

dr * = dr + (q — qo — ft ■ dr)n , 

where n is the normal on the wall that passes 
through r, qo is the distance of the particle 
at initial position r from the wall, and q is 
generated from pu(q, qo, dt, v, D q ) with v = n- 
F/rj and D q = D (compare relation ([8]) and 
Appendix |A| . 

We note that dt has to be adjusted such that A 3> C(|dr|) 
holds true with A being the smallest radius of curvature 
of the wall. 



B. Justification and discussion 



IV. HARD WALLS 

Our novel procedure to incorporate hard-body interac- 
tions into a Brownian dynamics simulation is now devel- 
oped for space-fixed walls. The walls might be curved so 
that their normal vectors attain a position dependence. 
The discussion will be carried out for point like Brownian 
particles. For extended spherical particles the distance to 
the hard surface is reduced by the radius of the particle. 

To develop the method independently of the used algo- 
rithm for the numerical integration we discuss the proce- 
dure on the level of the Langevin equation ([2]) formulated 
in terms of differentials with the tacit understanding that 
dr has to be seen as the numerically obtained displace- 
ment for the (discrete) time step dt. We also assume 



For the time step that has produced the unphysical 
configuration the displacement can be decomposed into a 
displacement perpendicular to the wall (along the normal 

ft), 

dfl = (ft ■ dr)n = —dt + V2DdW± , (25) 
V 

and a displacement parallel to the wall, 

dfu = df- (n-dr)n = - J -dt + V2DdW\\ . (26) 
?y 

The perpendicular component of the deterministic force 
is given by F± = [ft ■ F)n whereas the parallel compo- 
nent is -Fj| = F — F±. The components of the fluctuat- 
ing force are obtained correspondingly and are given by 



9 



dW± = (dW -n)n and dW\\ = dW~ (dW-ft)ft. These two 
components are uncorrelated due to the orthonormality 
of the perpendicular and parallel directions: 

(dWjMWjf ) = (dl^W(d^ - dWn v n p )) (27) 

= dtn a n fi {l~n ti n ti ) = 0. (28) 

The correlations need to be considered only for the cur- 
rent time step where an unphysical position is generated. 
The normal ft sets up a locally fixed coordinate system. 
However, one is not interested in its dynamical proper- 
ties so that its (possible) position dependence has not 
to be considered in a dynamical (co-moving) sense. For 
this particular time step the random forces in the two 
equations are those that appear at time t, for the cor- 
responding diffusion constants one has D± = D» = D. 
It is essential that the perpendicular and parallel com- 
ponents of the displacement decouple and decorrelate for 
this local coordinate, so that they can be considered in- 
dependently. 

Let us assume for the moment that we have an appro- 
priately small time step so that the typical displacement 
dr is small compared to the (smallest) radius of curvature 
of the point on the wall that corresponds to the normal 
ft£& Hence, the wall appears to be flat on the scale of dr. 
In this approximation, the perpendicular motion along 
the normal ft thus corresponds to a motion on an infinite 
half-line, whereas the motion along the parallel direction 
is unrestricted. The overlap with the wall then solely re- 
sults from the (one-dimensional) displacement along ft. 
A valid physical end position along that direction, i. e. 
the distance q from the wall after elapsed time dt, is now 
directly generated from <JSj> , with go being the separa- 
tion of the initial position r from the wall, D q = D and 

v = ft ■ F/rj. 

The parallel motion is uncorrelated with the perpen- 
dicular motion for the considered time step on the level 
of the fluctuating forces and therefore is unaffected by 
that modification. We note that the deterministic force 
is fixed for the considered time step and that the modi- 
fication accounts for the altered statistical properties of 
the random forces in the vicinity of a hard wall. Hence, 
the already proposed displacement dfu can be retained 
to avoid unnecessary computational efforts. The total 
modified physical displacement is then given by dr * 
(it ft) + dr|| with dq — q — go- Plugging in the originally 
proposed dr on gets the modified dr * = df+(dq—ft-dr)ft 
already given above. 

We note that the maximum allowable size of the in- 
tegration time step dt is limited by the curvature of the 
hard wall. However, the length scales of physical sys- 
tems for which Brownian dynamics is applied is typically 
mesoscopic so that walls are only weakly curved on the 
scale of typical displacements and thus the corresponding 
radius of curvature is large. To be more explicit, con- 
sider protein-protein association as a paradigmatic ex- 
ample. Here the extension of (space-fixed) obstacles is 
usually comparable to (or even larger than) the diame- 
ter of the proteins themselves, say 40 A. The range of the 



involved non-covalent forces is typically of the order of 
2 to 4A£i2ii2£ Therefore, the time step that has to be 
used to capture the variations in the force fields is al- 
ready small enough to account for the curvature of the 
obstacles as well. In summary, the local flatness of the 
boundary is necessary for our algorithm to be applicable, 
but this approximation can be controlled by the choice 
of the size of the time step very much in the same way 
as space-dependent variations of the forces are captured 
by sufficiently small time steps. 

Finally, we want to point out that in the force-free 
case v = 0, the solution ^ simplifies to p(q, t;qo) = 
Pi {q, t;q )+ pi (q, t;-q ) (since p 2 (q, t; q ) = pi (q, t; -q ) 
and P3(q,t;qo) = 0, compare consistent with what 
one would obtain from the so-called image method (see 
also Appendix [B]) . As a consequence, our new algorithm, 
which is constructed from expressions ([6]) and ([7]), be- 
comes equivalent (in a strict mathematical sense) to nu- 
merical schemes which replace invalid positions by an 
image end-point whenever a wall is encountered and ex- 
ternally applied forces are absent. A prominent example 
is the algorithm developed by Scala et al^ 

V. HARD SPHERES 

In the previous section physical systems with space- 
fixed hard walls have been considered. In many systems, 
however, movable objects with hard cores are present as 
well. In this section we describe, how the previously dis- 
cussed method can be extended to encounters of two hard 
spheres, which both move under the influence of deter- 
ministic forces and diffusion. 



A. Description of algorithm 

Consider two hard spheres with (possibly different) 
radii a,, friction coefficients r\i and diffusion constants 
Di = k B T/r]i whose dynamics is described by the two 
Langevin equations 

dri = vxdt + ^2LhdW! (29a) 

and 

dr 2 = v 2 dt + v / 2^dH> 2 . (29b) 

Here again the drift terms are denoted in terms of veloc- 
ities Hi = Fi(fi,t)/r]i for the two particles. 

Suppose now that the two spheres are at f\ and r 2 at 
time t and that we have at hand a standard integration 
scheme to update the particle positions for a time step 
dt. Then, the algorithm for generating valid displace- 
ments dr* and dr 2 * consistent with hard- wall interactions 
between these two particles is as follows: 

1. Use the standard integration scheme to generate 
'preliminary' propagation steps dr*i , dr2 . 
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2. Check whether or not the suggested end positions 
f\ + dfi and f 2 + dr 2 represent a physically valid 
configuration, i. e. whether or not the two spheres 
at positions fi + dfi and r 2 + dr 2 would overlap. 

a. If the two spheres do not overlap, the new 
configuration is accepted, drj" = dfi and 
df* = dfi. 

b. If the two spheres do overlap, valid displace- 
ments dfi*, df 2 * are constructed according to 
the following procedure: 

• Calculate the separation qq = |fa — fi| — 
(ai +02) of the surfaces of the two spheres 
and determine the normalized connection 
vector e= (f*2 — fi)/|f*2 — fi|- 

• Generate a position q from 
pn(q,qo,dt,v, D q ) with diffusion con- 
stant D q = D\ + D 2 and velocity 
v = (v 2 — fi) ' e (compare relation (|8|) 
and Appendix [AJ . 

• Use dq = q — q® to evaluate the modified 
physical displacements 

dfi* = dfi H — — ((dfa — dfi) • e — dq) e (30a) 

vi + m 

and 

d?^,* = d?"2 — — ((df*2 — dfi) ■ e — dq) e. (30b) 

vi + m 

The time step dt has to by chosen small enough so that 
ai + a 2 ^S> 0(|df2 — dfi |) is satisfied. 

B. Justification and discussion 

To generate a physical configuration without overlap 
the motion of the two Brownian particles is decomposed 
into a common center-of-friction motion and a relative 
motion of the two spheres. The center-of-friction R = 
(^lfi+^fi)/'?, T) = 771+772, obeys the Langevin equation 

dR = - (?7ifi + 7721/2) 
77 

+ - (mV2D~idW 1 + i l2 V^hdW 2 ) (31) 

= Vdt + y/2D R dW R (32) 
and the relative motion of r = f 2 — fi is governed by 

df* = (v 2 - vi)dt + ^/2D~ 2 dW 2 - y/2LhdWi (33) 
= vdt + y/2D^dW r . (34) 

Here, the abbreviations v := tj 2 — v\ and V 
^ (Vi^i + ^f^) = (Fi + F2) have been introduced. The 
fluctuation terms in (f3"2")l and (|34p can again be writ- 
ten in terms of un-correlated white noise terms dW R 
and dW r with associated diffusion coefficients Dr ~ 



k^T/i] = D 1 D 2 /(D 1 + D 2 ) and D r = D x + D 2 , respec- 
tively. The statistical properties of the fluctuation terms 
in (f3"2"|) and (|3^|) are readily obtained from the observation 
that the sum and the difference of two Wiener processes 
are again Wiener processes (as theses processes are Gaus- 
sian stochastic processes with vanishing mean). Due to 
the appropriate weighting of the contributions of the two 
particles in the definition of the center-of-friction, d\V R 
and dW r are uncorrelated, as can be seen directly by 
equating the cross correlations. 

As a consequence, any corrections in the relative mo- 
tion do not affect the propagation in center-of-friction R 
(and vice versa). The overlap of the two particles is the 
result of the unphysical propagation in the relative coor- 
dinate r so that this very motion has to be modified due 
to the hard-wall interaction, while the initially proposed 
displacement dR of the center-of-friction can be kept to 
avoid unnecessary computational efforts. This is similar 
to the situation of a particle in the vicinity of a wall where 
the normal and the parallel displacements decorrelate. 

The relative motion itself corresponds to the motion 
of a point-particle close to a space-fixed sphere of radius 
ai + a 2 . The treatment of the relative motion therefore 
follows the lines of reasoning already developed in Sec. 
IIV1 The relative displacement df can be further decom- 
posed into a (one-dimensional) displacement dfi along 
the connection line defined by the unit vector e := f/|fj 
and a displacement df|| in the plane perpendicular to e. 
Thus one obtains the equations 

dr ± = v ■ edt + y / 2D r dW r ■ e (35a) 

and 

df,, = (f - ujjdt + v/2Z^ (dW r - (dW r ■ e)ej (35b) 

with fj_ = (e • v)e. (Note that e is the normal of the 
effective hard surface that restricts the motion of the 
relative coordinate. The nomenclature of perpendicu- 
lar and parallel displacement has to be seen with respect 
to this spherical surface). The crucial point is that the 
perpendicular fluctuating term dW± — dW r ■ e and the 
parallel random term dl^i| = dW r — (dW r ■ e)e are again 
un-correlated (for the current time step where the over- 
lap is detected). The associated diffusion constants are 
D± = Dn = D r . The vector e is fixed for the considered 
time step and again defines a local coordinate system 
suitable for the decomposition. 

Assuming that the effective hard wall (that corre- 
sponds to a sphere of radius a\ + a 2 ) can be well ap- 
proximated as flat (which is, for example, well fulfilled 
in the typical case in point of protein-protein associa- 
tion, see Sec. HVj) . the relative motion along the direc- 
tion e then corresponds to the motion on an infinite half- 
line and the parallel motion is unrestricted in an infi- 
nite plane. Therefore, the described algorithm of Sec. [TT] 
for the one-dimensional system can be applied directly 
to generate a new physical position along e. The prop- 
agation in the plane perpendicular to e is not affected 
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by that modification and thus the corresponding initially 
obtained displacement can be kept. To satisfy the as- 
sumption that the hard surface is virtually flat the time 
step has to be chosen sufficiently small so that one has 
ai +a 2 -»0(\dr\)&. 

The new configuration without overlap of the spheres 
is obtained from the modified relative displacement by 



dr * = dry + dqe 



(36) 



with dfu = dr — (dr • e)e being the originally proposed 
parallel displacement (dr denotes the initially proposed 
relative displacement dr*2 — drj of the two spheres). Here, 
dq = q — q is fixed by q that is generated from distribu- 
tion (JSJ) with diffusion constant D q = D± = D\ + D2, 
velocity v = v ■ e — {V2 — V\) • e, initial separation 
qa = \T2 — r\ I — (ai +02) of the two spheres. Once dr * is 
determined the new absolute particle positions without 
overlap are calculated from the new displacements as 



and 



dr * = dR 



dr* = dR 



>)2 



m + m 



-dr 



m + m 



-dr*. 



(37a) 



(37b) 



Plugging in the explicit expressions for the various quan- 
tities one ends up with equations (|30aP and (j30b[) . 



VI. EXEMPLARY MODEL SYSTEMS 

In this section we numerically investigate the appli- 
cability and quality of the new algorithm in compari- 
son with the heuristic methods of Sec. IIIII We do not 
aim at a study of particular physical systems or phe- 
nomena, we rather intend to present exemplary problems 
that typically have to be dealt with in systems modelled 
within Brownian dynamics. To set up an explicit frame 
we look at two research areas from soft matter, namely 
the motion of colloids in flow fields and the interaction 
of proteins in solution. For simplicity we consider only 
two-dimensional systems. All numerical results presented 
below are obtained by using the Euler algorithm as stan- 
dard integrator. 



A. Colloids in flow fields 

The first characteristic system we look at is given by 
colloidal particles in stationary flow fields u(r). To begin 
with we consider a spherical particle of radius a = 1 /im 
and with diffusion constant D = 0.21/xm 2 /s (Stokes- 
Einstein diffusion at room temperature k B T room ss 
4 fN/im) suspended in the fluid flow around a space-fixed 
circular post in two space dimensions. This model system 
is taken to mimic a typical situation in microfluidic de- 
vices where particles encounter posts or walls in the chan- 
nel. Let us consider a flow field with a typical velocity of 



uq =10/im/s (typical velocities in microfluidic channels 
range from 1 /im/s to 1 cm/s) around a circular post of 
radius R = 10/itm at the origin of the coordinate system. 
For a time step dt = 0.01 s the typical diffusive displace- 
ment yj (r 2 ) — (r) 2 of a particle with radius a — 1 /im is 
0.1 /im. Thus we deduce that for a radius of R = 10 /im of 
the post the curvature of the hard surface should not limit 
the applicability of our method. Assuming no-slip bound- 
ary conditions, the velocity field of the two-dimensional 
stationary Stokes flow u{r) = (u x (x,y),u y (x,y)) around 
the post is given by 



u - 



y 



1 



R 2 



2uq In 



R 



and 



-uq- 



■ry 



1 



R 2 



(38a) 



(38b) 



in Cartesian coordinates (for a discussion of the Stokes 
approximation for small Reynolds numbers and its valid- 
ity in two space dimensions, see for example Refs. l37land 




x [(Jin] 

FIG. 5. Contour plot of the distribution of a spherical colloid 
(radius a = 1 /im, diffusion constant D = 0.21 /xm 2 /s) flowing 
around a circular obstacle of radius R — 10 /im. The par- 
ticle starts initially at r(0) = ( — 11.5 /im, /im) and evolves 
for t = 10 s. The Stokesian flow (|38[) is characterized by 
uo — 10 /im/s. The solid line represents the solution of the as- 
sociated Fokker-Planck equation. In the upper panel symbols 
display the numerical solutions obtained from integrating the 
Langevin equation with v(f) = u(f) for a time step dt = 0.01 s 
with circles corresponding to our algorithm and open (red) 
squares to the event-driven scheme. The lower panel shows 
the data from the rejection scheme for dt — 0.01 s (blue dia- 
monds) and dt = 0.001s (open green diamonds). The black 
region represents the obstacle, the grey shaded region is inac- 
cessible to the particle coordinate. 



Figure [5] compares the numerically determined distri- 
bution function of a colloid in flow field (13"81) with the 
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solution of the associated Fokker-Planck equation. The 
evolution of the particle as exhibited by the rejection 
scheme lags behind the expected behavior. To analyze 
this further we look at the mean first passage time the 
colloid needs to cross a certain threshold after passing 
the post. Figure [5] displays the corresponding results as 
a function of the integration step dt. It shows also the 
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FIG. 6. Left panel: Mean first passage time r a colloid of 
radius a = 1 /im needs to flow past a post of radius R = 
10 fim to the threshold at x = b = 11.5 (xm as a function 
of the integration time step. The colloid starts initially at 
r(0) = ( — 11.5 j-im, /tm). Our method is displayed by (black) 
circles, the event-driven scheme by (red) open squares and 
the rejection scheme by (blue) diamonds. For the numerical 
determination of the mean first passage time, see e. g. Refs. 
and |36| . Note that the collision frequency with the post 
drops from 3% to 1% for the shown range of dt. Right panel 
and inset: Tangential diffusion fa — (fdr • euj )/(2dt) during 
the first time step of size dt = 0.01s as a function of the 
initial distance qo = 1^(0) — (R + a)\ between the surfaces of 
the colloid and the post. The different schemes are displayed 
using the symbol code of the left panel. The expected value 
is 0.21 fim 2 /s. 



particle diffusion &2 = ( (dr • e[\ ) ) / (2dt) tangential to the 
post during the first time step as a function of the initial 
separation from the post's surface (eii denotes the tangen- 
tial unit vector to the surface). The tangential motion 
is not restricted by the presence of the hard surface so 
that ki is expected to be given by the diffusion constant 
D of the particle. For the rejection scheme, however, the 
tangential diffusion close to the wall is clearly reduced. 
This has the effect that the time the particle needs to flow 
around the obstacle is increased when measured from the 
rejection scheme. As can be seen in Fig. [6]this effect is re- 
duced for decreasing time steps as the collision frequency 
with the obstacle goes down as well. 

As a second typical case study we investigate the evo- 
lution of two colloidal particles in a prescribed shear flow, 
The velocity field at r = (x, y) is given by 



u(r) = aye x 



(39) 



where a is the shear rate and e x denotes the unit vector 
in ^-direction. Figure [7] depicts the distribution function 
of the relative coordinate of the two colloids for different 



elapsed times. Again numerical data is compared to the 
solution of the Fokker-Planck equation. 
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FIG. 7. Contour plot of the transition probability of the 
relative vector r — (x,y) — ?2~t\ of two colloids in a 
shear flow l|39[) with shear rate a — 10/s. The two particles 
with radii a\ = 1 /im, a-z — 0.5 fim and friction coefficients 
rn = 18.9 fNs/Atm, r] 2 = 9.45 fNs//xm (i.e. Di = 0.21 ^m 2 /s, 
D2 = 0.42/im 2 /s at room temperature) initially start at 
fi(0) = (0/im,0/im) and r*2(0) = (— 2 pim, 0.75 /mi), and 
evolve for t = 0.5 s (panels (a) and (b)). The solution of 
the associated Fokker-Planck equations are depicted by solid 
lines, the numerical results from integrating the Langevin 
equations with dt = 0.002 s are shown by symbols: (Black) 
circles represent our algorithm, open (red) squares correspond 
to the event-driven scheme (panel (a)), and (blue) diamonds 
to the rejection approach (panel (b)). Panel (c) shows the 
evolution of one contour line with time (0.5, 0.8, 1.0 s) in 
comparison to the results of the rejection scheme. The grey 
circle represents the region which is inaccessible to the relative 
coordinate of the two colloids. 



B. Protein-protein interaction 

A (fairly small) protein has a typical radius of approx- 
imately a = 2 ran leading to a diffusion constant D = 
0.1nm 2 /ns at room temperature k^T TOom w 4pNnm. 
Therefore, the typical time a diffusing protein needs to 
cross its size is of the order of 10 ns. The forces between 
two proteins are governed by non-covalent interactions 
which have a typical energy scale of 2 to 4 k B T mom and 
a range of 0.2 to 0.4 nm. Therefore, the velocity in the 
Langevin equation ([2]) is of the order of 5 nm/ns. A time 
st ep of dt = .001 ns has a typical diffusive displacement 
yj {f 2 ) — (r) 2 of approximately 0.02 nm and should be 
small enough to capture typical variations in the force 
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fields (0(|?7di|) ~ 0.005 nm) and to resolve the details of 
trajectories. 

To have an explicit example at hand we will look at the 
interaction of proteins such as lysozyme, which has a ra- 
dius of 1.8 nm and a diffusion constant of 0.12nm 2 /ns in 
aqueous solution. The force between two proteins stems 
from repulsive (screened) electrostatic and from attrac- 
tive short-ranged van-der Waals interactions. Modelling 
the lysozyme solution as a colloidal syste m 39 i 40 within 
the Derjaguin-Landau-Verwey-Overbeek theory the elec- 
trostatic interaction between two proteins of separation 
r = | ?2 — fi | is represented by the Debye-Hiickel contri- 
bution 

0dhM = — cxp (-«(r - a)) (40) 
r 

and the attractive van-der Waals part by 

A fa 2 a 2 (r 2 -o 2 \\ 

(41) 

The total interaction potential is then given by 

6(r) = [ </,DIi(r) + ^ HA(r) ' % r ~ ° + \ (42) 
^ w 1 oo, if r < a + 6 y ' 

from which one obtains the force Fi = —Vi4>(r) on par- 
ticle i = 1,2. Here k denotes the screening length, a 
is the radius of the effective particle (i. e. the sum of the 
two radii) and S is the Stern layer thickness so that the 
excluded hard-core radius is given by Rue = u + 5. 

For two proteins that experience only their mutual in- 
teractions, the motion of the common center-of-friction 
completely decouples from their relative motion. Let 
us first look a this relative motion, which then corre- 
sponds to a driven diffusive motion of a point-particle 
around a sphere. The transition probability for the rela- 
tive coordinate in the vicinity of the excluded hard-core 
region is shown in Fig. [5] The time step is chosen to 

be dt = 0.001ns so that O (V(r 2 ) - (r) 2 ) - 0.02 nm. 

Therefore, we expect that the curvature of the excluded 
hard-core sphere does not lead to considerable system- 
atic deviations from the exact transition probability. The 
numerically obtained transition probabilities for one in- 
tegration step are compared with the 'exact' transition 
probability that has been calculated by solving the cor- 
responding Fokker-Planck equation. 

C. Mean first passage time 

As an example of a derived quantity we look more 
closely at the mean first passage time for two different 
driven diffusive processes where hard-wall interactions 
play an important role (see also Fig. [S] for a first brief 
example). Although these cases capture typical situa- 
tions one encounters in simulations of soft matter sys- 
tems, they are simple enough so that the mean first pas- 
sage time obtained from the numerical integration of the 
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FIG. 8. Contour plot of the transition probability of the rel- 
ative vector of two proteins. Its evolution corresponds to 
a point particle with D = D\ + D2 = 0.24nm 2 /ns. The 
relative coordinate is initially r(0) = (4.025 nm, nm) and 
evolves for the time t = 0.001 ns. The force (for the rel- 
ative component, compare relation ()34p ) is obtained from 
(|42|) with K = 4 fceTroom nm, k = 0.02 nm" 1 , a = 3.6 nm, 
A — 3fc B Troom an d 5 = 0.4 nm (compare e.g. Ref. 1401 ) mak- 
ing up an excluded hard-core region of radius -Rhc = 4nm. 
The numerical results are obtained from one integration step 
(i.e. dt = 0.001 ns) for which the encounter frequency is ap- 
proximately 12.7%. Exact curves obtained by solving the 
Fokker-Planck equation are displayed by solid lines. Left 
panel: (black) filled circles for our new algorithm, open (red) 
squares for event-driven scheme, for symmetry reasons only 
positive y axe shown. For the inner contour levels results from 
our method are shown separately in the upper right panel, 
those from the event-driven are shown in the middle panel on 
the right. The lower right panel exhibits the same data for 
the rejection scheme (blue filled diamonds). The grey shad- 
ing represents the region which is inaccessible to the relative 
coordinate. 



equations of motion can be compared with an analytical 
solution. The latter is obtained from the mean escape 
time from an interval [a, b] with a reflecting boundary at 
a and an adsorbing boundary at b, 

b q 

r(«Zb) = ^ / dqcMPHq)) J d</cxp(-/3<%')) (43) 

90 a 

for a particle at inverse temperature j3 = l/k B T moving 
in a potential $((7), and starting at q & 

In the first situation we continue the example from 
the previous section and consider the mean time the two 
proteins need to depart a certain distance b. This is re- 
lated to a thermally activated process where the reaction 
coordinate is the separation r = \r\ = |f^j — r*i|. Tak- 
ing the interaction potential (|4"2")l the problem is radially 
symmetric suggesting a transformation to polar coordi- 
nates. Due to the isotropy of the fluctuating forces in 
the Langevin equation the corresponding random forces 
in orthogonal polar coordinates are uncorrelated. There- 
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fore, the escape time is only determined by the evolution 
equation for the separation r which can be obtained from 
by applying Ito's formula^: 

dr = —l-l$(r)d* + ^2D~ T dW r (44) 
rj r dr 

where r) r = k B T/D r is the friction associated with the 
relative coordinate (D r = D\ + D2) and dW r is an incre- 
mental Wiener process. The effective potential is given 
by $(r) = 4>(r) — fc B Tlnr in two space dimensions. Fig- 
ure [9] depicts the exact values of the escape time as a 
function of the threshold b together with the numerically 
obtained values. Before the particle crosses the threshold 
at b many integration steps have to be carried out. In 
addition, a small time step dt has to be chosen to capture 
the variations in the force field so that only a small frac- 
tion of encounter events will occur. It is thus expected 
that the heuristic algorithms will produce good results. 
Figure [9] indeed confirms this expectation. This example 
points out that the rejection scheme may produce accu- 
rate results for certain observables. However, this does 
not necessarily mean that the rejection procedure is reli- 
ably applicable in general situations, as already demon- 
strated with the examples above, and as will be under- 
lined by the following findings for the mean first passage 
time of a particle close to a wall. 
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FIG. 9. Escape time r from the potential (|42p as a function 
of the location b of the (radial) threshold. The same potential 
parameters as in Fig. [8]are used with r(0) = (4.025 nm, nm) 
and dt = 0.001 ns. With increasing b the percentage of 
time steps that involve a 'collision' with the hard-core region 
drops from approximately 13% to 2% for the shown data. 
The numerical data is shown by (black) circles for our pro- 
posed method, by (red) open squares for the event-driven 
scheme and by (blue) diamonds for the rejection scheme, 
the solid curve represents the exact escape time evaluated 
from (|43[) with a — i?HC = 4nm, go = 4.025 nm and 
<£(?) = <f)(q) — fe B Tlnr (see (142 1) and main text). The mean 
escape time was measured from the simulations as described 
in Refs. and I3H 

As a second example we consider a (colloidal) particle 
which hits a hard wall (for instance, a wall or obstacle 



in a microfluidic device or an organelle in a cell etc.). 
The particle is driven by a deterministic force with a 
non-vanishing component perpendicular to the wall (for 
example due to an electric or magnetic field). For sim- 
plicity, we restrict ourselves again to two dimensions and 
assume that the wall located at x = is infinitely ex- 
tended along the y direction such that it confines the 
particle to x > and y £ M. The particle is subject 
to the force F = (— /, /). In Fig. [TU] we show the mean 
time a particle of radius 1 /jm starting at r (0) = (0, 2 fim) 
needs to 'slide along the wall' by a distance of 10 parti- 
cle radii until it crosses the line at b = y = 10 /iin. The 
results obtained for the different integration schemes are 
compared to the exact solution calculated from (|4"5|) with 
Q(q) = —fq to be r = brj/f. While our new integra- 
tion scheme yields the exact passage times, the rejection 
scheme dramatically fails for too large / as it artificially 
slows down the particle motion due to prevailing rejec- 
tion steps until it practically stops completely. This goes 
in hand with a considerable increased computing time, 
see inset of Fig. 1101 The event-driven scheme produces 
very good results. However, it requires (exponentially) 
increasing computation time in the case of high collision 
frequency and strong forces into the wall as many realiza- 
tions of the random force for the propagation step away 
from the wall will produce invalid end positions. 



D. Computational efficiency 

To analyze the computational efficiency of our algo- 
rithm in comparison to the rejection and event driven 
scheme in some more detail and to get an impression 
of its performance in realistic systems, we come back to 
the above example of the two interacting proteins. As in 
realistic applications, we simulate the evolution of the in- 
dividual coordinates of the two interacting particles, i.e. 
we do not decouple the center-of-friction motion. The 
event-driven scheme requires the determination of the 
collision time of the two spheres once an overlap has been 
detected. This collision time can be determined analyt- 
ically for this special case by solving a quartic equation 
for non-vanishing forces (for zero force it reduces to a 
quadratic equation). Having the possible generalization 
of an arbitrarily shaped rigid particle in mind, however, 
we calculated the collision time by an iterative scheme 
which was used by Tao et alJ^ They determined the col- 
lision time by subdividing the integration step, which led 
to a particle overlap, into 100 sub-steps to successively 
propagate the particles along the "Euler path" until the 
collision is detected. Alternatively , the collision time can 
be determined by using nested intervals. Our method 
also requires the (numerical) evaluation of a root when 
determining the modified normal component of the rela- 
tive motion (see Appendix [A| . However, the motion per- 
pendicular to the normal component and, in particular, 
the evolution of the center coordinate need not be altered 
and are indeed kept unchanged by the algorithm, so that 
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FIG. 10. Mean first passage time of a particle (radius 1 /im) 
close to a wall (located at x = 0) as a function of v = f/r). 
It moves from r(0) = (0, 2 (im) to the line y = 6 = 10 ^m 
driven by the force F = (— /, /) (see main text). (For fixed / 
the convergence as a function of dt is similar to Fig. [U not 
shown.) The numerical data is shown by the (black) line for 
our method, which is indistinguishable from the exact passage 
time t = b/v evaluated from (|43[) with a — — oo, b = 10 /im, 
go = and 3?(g) = — /g = —r/vq. The (red) dashed line 
represents the event-driven scheme and the (blue) dotted line 
the rejection scheme (time step dt = 0.01 s in all cases). The 
mean escape time was measured from the simulations as de- 
scribed in Refs. and HH. The inset shows the computation 
time consumed by the different algorithms. (See Sec. IVIDI 
on how the collision time is determined for the event-driven 
scheme. Upper dashed curve: stepwise propagation; lower 
dashed curve: nested intervals.) For larger driving force, the 
computation times of the rejection and event-driven scheme is 
an order of magnitude longer as compared to our algorithm. 
For the shown data the 'collision' frequency increases from a 
few per cent to around 90%. 



only a one-dimensional motion is affected. We therefore 
expect our algorithm to be faster than the event-driven 
scheme where all components are simultaneously affected 
when accounting for hard bodies. For the typical param- 
eters of the exemplary protein system from above, see for 
example Fig. [8l we indeed found that the required com- 
puting time of our method is smaller than the time for the 
event-driven scheme. Not surprisingly, the discrepancy 
increases for increasing encounter frequencies w. For the 
data presented in this work we typically observe a differ- 
ence of 10 % to 30% in computational speed. The rejec- 
tion scheme, on the other hand, is typically 20% faster at 
identical time steps. However, to diminish the systematic 
errors due to encounter events with hard walls the time 
step has to be reduced accordingly by at least one magni- 
tude, thus requiring considerably increased computation 
time. 

For a more quantitative analysis we investigate the 
necessary computational efforts to achieve a prescribed 
accuracy. To this end we again use the passage time 
problem of a colloid gliding along a wall of the previous 
subsection as for this example analytical results can be 
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FIG. 11. First passage time distribution for the system de- 
scribed in Fig. [10] with v = 1 /im/s. Upper left panel: log- 
arithm of the simulated distribution of our method (black 
circles) and of the event driven scheme (open red squares) to- 
gether with the exact curve (solid brown line) for dt = 0.1 s; 
upper right panel: same comparison for the rejection scheme 
(blue diamonds) and in addition curves for smaller dt (thin 
black lines, dt = 10" 4/3 , 10~ 5/3 , 10~ 2 , 10 _7/3 s). The middle 
panels show the relative deviation of the simulated mean first 
passage time T s ; m from the exact value r = 10 s, the lower 
panels show the Jensen-Shannon divergence (|46[) . Both quan- 
tities are shown in double-logarithmic depiction as a function 
of the time step (left hand side) and the computation time 
(right hand side, the number of realizations is fixed to 10 7 , 
dt is varied, nested intervals are used for the event-driven 
method). Note the different order of magnitude for the differ- 
ent methods (our method: black circles; event driven: open 
red squares; rejection: blue diamonds). The horizontal brown 
line indicates a prescribed accuracy of (r s i m —t)/t = 1 %. 
For this setting the event driven-scheme exhibits a rather large 
encounter frequency of about 50 % with the wall and needs 
by a factor three more computation time than our meth- 
ods. For dt = 0.01 s (compare Fig. I10[) an accuracy of 
( T sim — r )/ T ~ 0.5 % is achieved, the event-driven scheme 
then produces an encounter frequency of 11 % and needs 
roughly 55% more computation time. In contrast, the re- 
jection scheme shows a very slow 'convergence'. 



obtained. For the system described in Fig. QJJ] (starting 
point yo = 0) the first passage time distribution Q(t) can 
be calculated^ to be 



Q(t) 



\b\ 



V^nDt 3 



exp 



(6 - vtf 
ADt 



(45) 



In addition to analyzing the relative deviation of the sim- 
ulated mean first passage time r S i m from the analytical 
result t = brj/ f we will compare this exact distribution Q 
with the simulated curve Q S i m by determining the Jensen- 
Shannon divergence^ as a measure for the distance be- 
tween Q and Q s i m . The Jensen-Shannon divergence is 
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defined by 

oo 

J =2l^/ dt ( QWb S| +0 - (t)ln %f) 



(46) 

where M = (Q + Q s i m )/2. J can take on values between 
and 1 with indicating identical distributions^ The 
findings for the system of Fig. [TU] are shown in Fig. QT] 
Surprisingly the event-driven scheme produces results of 
the same accuracy as our new methods for the same time 
step. Therefore, the computational difference is solely de- 
termined by the effort for a single integration step which 
amounts to typically 10-30 % for moderate encounter fre- 
quencies with the wall smaller than 10 %. In accordance 
with our general observation the data in Fig. [11] point 
out again that a rather small time step has to be cho- 
sen for the rejection scheme to match the accuracy of the 
other methods. 



VII. CONCLUSIONS AND OUTLOOK 

The explicit form of the transition probability for a 
Brownian particle in a constant force field on the one- 
dimensional half-line shows that multiple encounters with 
hard bodies may lead to contributions to that transition 
probability that cannot be generated by a direct numer- 
ical integration of the Langevin equations. This results 
in unsystematic and uncontrolled errors in heuristic in- 
tegration schemes, even for very small time steps. For 
an algorithm whose accuracy is to be systematically con- 
trolled by the size of the discretization time step, the 
recourse to the transition probability seems unavoidable. 

In this article we have developed a novel algorithm to 
account for hard-body interactions in Brownian dynam- 
ics simulations with controllable discretization error, and 
compared its performance with heuristic schemes from 
the existing literature. The central idea is that, once a 
hard-body interaction has been detected during the nu- 
merical integration of the Langevin equations of motion, 
the component of the particle displacement (s) involved 
in the 'collision' is generated directly from the solution 
for the transition probability in the presence of a hard re- 
flecting wall. We demonstrated that the aforementioned 
transition probability for the one-dimensional half-line 
can be used to construct an approximate solution for 
general many-dimensional physical systems with hard- 
body interactions. The algorithm decomposes the par- 
ticle motion into parts that are unaffected by the hard- 
body interaction and into an 'affected' part along which 
the hard-wall 'collision' occurs. The latter corresponds 
to the relative motion of the collision partners and is 
reduced a one-dimensional motion on a half-line with a 
hard, reflecting boundary at the origin (the point of 'colli- 
sion'). For this setting, the general time-dependent tran- 
sition probability is known analytically, and can be used 
to generate relative displacements which are consistent 



with the hard-body interactions and with the statistical 
properties of the underlying stochastic evolution. 

This principle can be extended to more general situa- 
tions, for example, where structured particles with a non- 
spherical shape are present. In general, the decomposi- 
tion of the full particle displacement (s) has to be adopted 
in such a way that the components 'unaffected' by the 
hard-body interaction and the 'affected' component sta- 
tistically decouple from each other (for the considered 
time step), so that the hard-body 'collision' can be cor- 
rected without interfering with the rest of the particle 
motions. The 'unaffected' displacement is thus retained 
(to save computation time) and together with the cor- 
rected 'collision' part is used to reconstruct the overall 
displacements of the involved particles. In this work we 
developed this procedure for the case of interacting spher- 
ical particles, the generalization to more complicatedly 
structured rigid particles will be presented in a separate 
article. Furthermore, situations where multiple particle 
overlaps occur, for example in dense systems, and where 
secondary overlaps due to the applied modifications show 
up will also be treated elsewhere. 
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Appendix A: Generation of random numers distributed 
according to |[8]| 

One of the central building blocks of the described al- 
gorithm is the generation of a random number q that 
is distributed according to <JSj> . This can be achieved 
by using a general transformation method (see e. g. § 
or I44I ). Consider a distribution function p and define 
F(q) = Jq dq'p(q'). If the number x is uniformly dis- 
tributed on [0, 1] then q = F~ 1 (x) is distributed accord- 
ing to the distribution p. Here F -1 is the inverse of F. 

For the distribution function p{q) = pn(q, qo, dt, v, D q ) 
with ph given in ([8]) the corresponding integral is 




Due to the monotony of F the equation x = F(q) has 
precisely one solution q for given x. As to the best of our 
knowledge there is no closed expression for the inverse 
of F, q is obtained by solving x = F(q) numerically us- 
ing, for example, Newton's methods, bisection methods 
or combinations thereof. In this work Brent's scheme^ 
was used from the GNU scientific library^. 
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Appendix B: Harmonic oscillator on the half-line 

Let us consider a Brownian particle on a half-line in a 
harmonic potential with spring constant k, i. e. the total 
potential is given by 

m = { h - kq2 ' '% q -l (bi) 

YyHJ I co, if q < v ; 

with k > 0. Introducing the characteristic time r = 77/fc 
the corresponding overdamped Langevin equation reads 

dq = --dt + sj2D~ q dW (B2) 
T 

with a reflective boundary at q = 0. For the unre- 
stricted harmonic potential on the real axis the transi- 
tion probability p(q, t; go) of the Ornstein-Uhlenbeck pro- 
cess (|B2|) is analytically known^ for initial condition 
p(q,t = 0;q Q ) = S(q - g ): 

(B3) 

where 

s = s(t) = cxp (~J ■ ( B4 ) 

This solution can be used to construct the transition 
probability for the harmonic oscillator (|B1[I and reflecting 
boundary condition on R>o by using the image method: 

P(q, t; qo) = Pov(q, t; q ) + Pov(-q, t\ go) • (B5) 

Figure [T^] compares this exact solution with the nu- 
merically obtained transition probability using the algo- 
rithms discussed in the main text. To account for the 
variations in the force field the integration time step is 
chosen to be dt = O.Olr. This time step is already small 
enough that the solution from the event driven scheme 
is of similar numerical accuracy as our new integration 
scheme. For the rejection method, however, there are 
still deviations close to the wall due to the insufficient 
treatment of the reflecting boundary condition. 

The transition probability of the unrestricted 
Ornstein-Uhlenbeck process on R is a Gaussian. 
Therefore, an exact integration algorithm for any dt 
exists^: 

q(t + dt) = q(t)s(dt) + yrAjU - [s(dt)] 2 )G (B6) 

with G being a random Gaussian number of mean zero 
and variance one and s(dt) given in (|B4|) . Due to the 
symmetric structure of the transition probability (|B5j) 
one can readily construct an exact update scheme to gen- 
erate valid positions q* for the potential (|B1[) with a hard 
wall at the center as well: 

1. Calculate a 'preliminary' position q(t + dt) accord- 
ing to (|B6|) by a standard integration algorithm. 




FIG. 12. Transition probability of a particle on a half-line 
in the harmonic potential (|B1[I for the parameters r = 1.0, 
D q — 1.0 and initial position qo — 1.0. In the left panel the an- 
alytical solution (compare (|B5[> and (|B3|) ) is displayed for two 
times, namely t — 0.1 (solid line) and t — 0.2 (dashed line). 28 
For the generation of the numerical data the rejection method 
and our scheme is used with a time step dt = O.Olr = 0.01. 
(The results from the event-driven algorithm are not displayed 
as they are indistinguishable from our method). The upper 
right panel shows a blow up for t = 0.1, the lower one for 
t — 0.2 (same magnification factor). For our new method 
(red dashed lines) there are still tiny deviations which are re- 
lated to the finite discretization size of the time step. The 
rejection method (blue dotted lines) suffers from errors close 
to the wall due to the insufficient representation of hard-wall 
interactions. 



2 a. If this proposed q is in the physical domain, i.e. 
q > 0, it is accepted as the new position, q* = q. 

2 b. If this proposed position is unphysical, i.e. q < 0, a 
new physical position q* > is obtained by simply 
reflecting g, i.e. g* = — g. 

We finally remark that this procedure, the construction 
of valid end positions by reflection, is possible because the 
exact solution can be obtained by the image method. For 
k the system reduces to free diffusion on a half-line 
and the described method reduces to the one discussed 
in Ref. [Tt] in the special case of one space dimension. 
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